Problem statement: In this project, we tried to make classifications of the trails based on their length and change in elevation. We also made one function that would simplify the process of classification. To generate this analysis, we used boundary, buildings, contours_3m, and trails layers from macleish package and contours250k layer from MassGIS.

#devtools::install_github("beanumber/macleish")
library(tidyverse)
library(sf)
library(macleish)
library(leaflet)
library(dplyr)
library(lwgeom)
library(ggplot2)
elevation <- st_read("contours250k/CONTOURS250K_ARC.shp")
plot(st_geometry(elevation))

# add data source

We first plotted the elevation map of Massachusetts illustrating the number of feet or meters the state rises above sea level. Since the graph plots contours on the map, the denser the contour lines, the higher the place above the sea level. The highest point in Massachusetts is in the northwestern part of the state, in Berkshire County, while the eastern third of Massachusetts is relatively lower than the other area.

elevation <- st_transform(elevation, 4326)
elevInter <- 
  elevation %>%
  st_intersection(macleish_layers[["boundary"]])

ggplot() + 
  geom_sf(data = elevInter, aes(color = factor(CONTOUR_FT), colors = "Set1")) +
  geom_sf(data = macleish_layers[["trails"]], aes(color = factor(name), alpha = 0.8, lwd = 1.2)) +
  theme(legend.text=element_text(size=15))

To rate the difficulty of the trail, we wanted to check the change of elevation. We made a contour map of the MacLeish Field Station, using different colors of lines to represent different contours and plotted the trails on the contour map. In this graph, each intersection of contour lines and trails represents a 30 ft change in elevation of that trail. To see intersections more clearly, we changed the alpha and the size of the lines so that pathways can stand out from the contour lines. We created a leaflet map to make an interactive visualization of elevation and trails as follows.

leaflet() %>%
  addTiles() %>%
  addPolylines(data = st_geometry(elevInter), weight = 2) %>%
  addPolylines(data = st_geometry(macleish_layers[["trails"]]), color = macleish_layers[["trails"]]$color) %>%
  addProviderTiles(providers$OpenStreetMap)
st_transform(elevation, 4326)
elevation %>%
  st_intersection(macleish_layers[["trails"]])
#change the color name in Macleish to lower case so that r can understand it
newTrails <- macleish_layers[["trails"]] %>%
  mutate(colorName = tolower(color))

markers <- 
  macleish_layers[["trails"]] %>%
  group_by(name) %>%
  filter(geometry == min(geometry))
#interactive map with trails in Macleish
leaflet() %>%
  addTiles() %>%
  addPolygons(data = macleish_layers[["trails"]], label = ~name) %>%
  addPolylines(data = newTrails, weight = 3, color = newTrails$colorName) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  #a topographic imagery map layer
  addProviderTiles("Esri.WorldTopoMap", group = "Topo") %>% 
  #a satellite imagery map layer
  addProviderTiles("Esri.WorldImagery", group = "Sat") %>%
  #add the control to switch between the two
  addLayersControl(baseGroups = c("OpenSreetMap", "Topo","Sat"),
                   options = layersControlOptions(collapsed = F))

We also wanted to look at the length of each path. To visualize this, we plotted trails on the leaflet map on a topographic imagery map layer, a satellite imagery map layer and an open street imagery layer, and added the control to switch between those three imagery layers.

#calculate each part of the length of the trail
length <- 
  macleish_layers[["trails"]] %>%
  group_by(name) %>%
  st_length()

#make a new column in newTrails data frame
newTrails <- 
  newTrails %>%
  mutate(length = 1)

#change the value in the length column to the corresponding value in length set
for (k in 1:nrow(newTrails)){
  newTrails$length[k] <- length[k]
}

#get the total length of each trail
newTrails <- 
  newTrails %>%
  group_by(name) %>% 
  summarise(length = sum(length)) %>%
  arrange(desc(length))
trailInter <- 
  elevation %>%
  st_intersection(macleish_layers[["trails"]])

trailInter <- 
  trailInter %>%
  group_by(name) %>%
  summarise(ele = ifelse(n_distinct(CONTOUR_FT) != 1, sd(CONTOUR_FT), 0)) %>%
  arrange(desc(ele))

We wanted to look at the exact number of the length and the change in elevation to get a more accurate rating. Therefore, we generated the table newTrails to calculate the total length of each path and the table trailInter to calculate the change in elevation. For the elevation change, we computed the standard deviation of the elevation of each trail. Higher standard deviation represents more changes in elevation and higher level of difficulty, while lower standard deviation represents less changes in elevation and lower level of difficulty. If a trail doesn’t show up in the trailInter table, we would regard it as there is no intersection between the contour lines and that trail which means there is barely no change in elevation.

joinTrail <- 
  macleish_layers[["trails"]] %>%
  st_join(trailInter, join = st_crosses)
## although coordinates are longitude/latitude, st_crosses assumes that they are planar
newTrailsNone <- 
  st_set_geometry(newTrails, NULL)
trailInterNone <- 
  st_set_geometry(trailInter, NULL)

joined <- 
  newTrailsNone %>%
  left_join(trailInterNone, by = "name")

joined$ele[is.na(joined$ele)] <- 0

joined <- 
  joined %>%
  mutate(difficulty = 0.4 * as.numeric(length) + 0.6 * as.numeric(ele)) %>%
  arrange(desc(difficulty))

joined$ID <- seq.int(nrow(joined))

joined <- 
  joined %>%
  mutate(level = ifelse(ID < 4, "hard", ifelse(ID > 6, "easy", "moderate")))
macleish_layers[["trails"]] <- 
  macleish_layers[["trails"]] %>%
  inner_join(joined, by = "name")
pall <- colorNumeric(palette="viridis",
                     domain= elevInter$CONTOUR_FT)
palTrail <- colorFactor(c("red", "blue", "black"), domain = macleish_layers[["trails"]]$level)

leaflet(data=macleish_layers[["trails"]]) %>%
  addTiles() %>%
  addPolylines(data=elevInter, 
               color = ~pall(CONTOUR_FT), 
               group = "Contours Lines", weight = 1) %>%
  addLegend("bottomleft", pal = pall, 
            values = ~elevInter$CONTOUR_FT, 
            title = "Macleish Contour Line in feet") %>%
  addPolylines(data = macleish_layers[["trails"]], weight = 2, label = ~name, color = ~palTrail(level), group = "Trails") %>%
  addLegend("bottomright", pal = palTrail, 
            values = ~macleish_layers[["trails"]]$level, 
            title = "Difficulty levels for biking") %>%
  addMeasure() %>%
  addControl("Macleish Trails Ranked by Difficulty Levels", position = "topleft") %>%
  addMiniMap(
    toggleDisplay = TRUE) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Topo") %>% 
  #and a satellite imagery map layer
  addProviderTiles("Esri.WorldImagery", group = "Sat") %>%
  #add the control to switch between the two
  addLayersControl(baseGroups = c("OpenSreetMap", "Topo","Sat"),
                   options = layersControlOptions(collapsed = F)) %>%
  addEasyButton(easyButton(
    icon="fa-crosshairs", title="Locate Me",
    onClick=JS("function(btn, map){ map.locate({setView: true}); }")))
pall2 <- colorNumeric(palette="viridis",
                     domain= macleish_layers[["contours_3m"]]$ELEV_M)
palTrail <- colorFactor(c("red", "blue", "black"), domain = macleish_layers[["trails"]]$level)

leaflet(data=macleish_layers[["trails"]]) %>%
  addTiles() %>%
  addPolylines(data=macleish_layers[["contours_3m"]], 
               color = ~pall2(ELEV_M), 
               group = "Contours Lines", weight = 1) %>%
  addLegend("bottomleft", pal = pall2, 
            values = ~macleish_layers[["contours_3m"]]$ELEV_M, 
            title = "Macleish Contour Line in meters") %>%
  addPolylines(data = macleish_layers[["trails"]], weight = 2, label = ~name, color = ~palTrail(level), group = "Trails") %>%
  addLegend("bottomright", pal = palTrail, 
            values = ~macleish_layers[["trails"]]$level, 
            title = "Difficulty levels for biking") %>%
  addMeasure() %>%
  addControl("Macleish Trails Ranked by Difficulty Levels", position = "topleft") %>%
  addMiniMap(
    toggleDisplay = TRUE) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Topo") %>% 
  #a satellite imagery map layer
  addProviderTiles("Esri.WorldImagery", group = "Sat") %>%
  #add the control to switch between the two
  addLayersControl(baseGroups = c("OpenSreetMap", "Topo","Sat"),
                   options = layersControlOptions(collapsed = F)) %>%
  addEasyButton(easyButton(
    icon="fa-crosshairs", title="Locate Me",
    onClick=JS("function(btn, map){ map.locate({setView: true}); }")))

We then joined newTrails and trailInter together, caluculating the difficulty following the formula \(diffculty = 0.4 \cdot length+ 0.6\cdot ele\) Considering it’s more diffcult to ride uphills than the flat, we assigned 0.4 for length and 0.6 for the standard deviation of elevation. We ranked the difficulty, for the three trails with the highest degrees of difficulty, which are Snowmobile Trail, Eastern Loop, and Western Loop, we rated them as in the level of “Difficult.” For the three trails with the lowest level of difficulty, which are entry trail, Driveway, Easy Out, we classified them in the level of “Easy.” And the rest of three trails, which are Poplar Hill Road, Porcupine Trail, and Vernal Pool Loop, we rated them in the level of “Moderate.”

  1. Github link: https://github.com/ChristinaLyu/miniProject3.git↩